Probability distribution of Majorana end-state energies in disordered wires 
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One-dimensional topological superconductors harbor Majorana bound states at their ends. For 
superconducting wires of finite length L, these Majorana states combine into fermionic excitations 
with an energy eo that is exponentially small in L. Weak disorder leaves the energy splitting 
exponentially small, but afi^ects its typical value and causes large sample-to-sample fluctuations. 
We show that the probability distribution of eo is log normal in the limit of large L, whereas 
the distribution of the lowest-lying bulk energy level £i has an algebraic tail at small ei. Our 
findings have implications for the speed at which a topological quantum computer can be operated. 
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Introduction. — Quantum bits based on topologically 
protected states promise a platform for error-free quan- 
tum computation [l}^^. Since information is stored in 
states with a topologically protected degeneracy, qubits 
that rely on this principle are immune to local external 
sources of decoherence. In practical realizations, how- 
ever, the energy splitting of the topological qubit is not 
exactly zero, because of finite-size effects. This poses a 
restriction on the speed at which a quantum computer 
must be operated: Operations have to be performed in 
a time that is short in comparison to the inverse energy 
splitting, but long in comparison to the inverse excitation 
gap for (non-topological) excitations. 

A particularly promising possibility to realize topo- 
logically protected zero-energy states is found in one- 
dimensional spinless p-wave superconductors, which are 
known to have zero-energy Majorana fermion states at 
their ends [HE]. Although Majorana excitations are in- 
sufficient to build a universal topological quantum com- 
puter, their implementation may considerably reduce the 
minimum required accuracy of qubit operations. There 
are several proposals for the experimental realization of 
such topological superconducting wires In some 

of these, one-dimensional wires can be brought into an 
alternation of topological and non-topological domains, 
with Majorana bound states at the domain boundaries 
[3 mm], while the location of the Majorana states can 
be controlled via gate voltages or magnetic fields [5] [13] . 

Experimental realizations necessarily involve topolog- 
ical domains of finite length L, as well as disorder. For 
finite L, the Majorana end states fuse into fermionic ex- 
citations at a finite energy Eq that is exponentially small 
in L/^, where ^ is the superconductor coherence length 
[5]. In disordered wires, this sets a lower bound for the 
speed of qubit operations which is sample specific. On 
the other hand, disorder is known to cause a Lifschitz 
tail of localized states below the gap A [51 [T3]. Since 
operations with Majorana states require that they are 
transported through the quantum wires [13], the lowest- 
lying bulk state of energy ei provides an upper bound for 
the speed of qubit operations. In view of possible exper- 



imental applications and their limitations, it is essential 
to know the full probability distribution of the energies 
Eo and El. This problem is addressed in this paper. 

Solitons in Dirac equation with random mass. — We 
first consider the Dirac Hamiltonian with random mass, 

H ^ VFpa^ - A{x)ax ~V{x)ax, (1) 

as a simple model for a topological superconductor with 
Majorana end states. Here vp is the Fermi velocity, A 
the effective superconducting gap, and and are 
Pauli matrices. The disorder potential V{x) is taken ac- 
cording to a Gaussian distribution with zero mean and 
correlator {V{x)V{x')) = jS{x — x'), corresponding to 
the mean free path / = vpT — Vp/j in the normal state. 
(Long-range correlated disorder which breaks the system 
into topological and nontopological regions has been con- 
sidered in Ref. [13 [HI) The Hamiltonian ([I]) arises as a 
low-energy effective Hamiltonian for semiconductor wires 
with strong spin-orbit coupling in proximity to a conven- 
tional superconductor and in the presence of a magnetic 
field [T7]. In that case, the effective gap A in Eq. ([T]) 
is the difference between the applied magnetic field and 
the proximity-induced superconducting gap in the ab- 
sence of a magnetic field. The same Hamiltonian arises 
in a number of other contexts, such as fermions on a lat- 
tice with random hopping amplitudes [TB], narrow-gap 
semiconductors [TH [50] , or organic molecules [Hj . What 
appears here as a pair of Majorana end states is referred 
to as a "soliton-anti-soliton pair" in these contexts [22] . 

We describe a topological domain of length L by set- 
ting A(a;) — — oo for x < and x > L and A(a::) — A for 
< X < L. Since the system is fully gapped for 2: < 
and X > L, any states contributing to the integrated 
density of states A''(e) = d£'v{E') must be localized in 
the topological domain at < a; < L. In the presence 
of disorder, N{e) is a random quantity with probability 
distribution Pg{N), which is related to the probability 
distribution Pj{£) of the energy level Ej {j = 0, 1, 2, . . .) 
through the equalities 

Pde) = -§-^i2P^(f)- (2) 
j'=o 
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Thus the distribution functions of the Majorana end- 
state energy eo and of the lowest bulk state en- 
ergy El obey po{£) — ~dPe{Q)/de and pi{e) — 
-d[P^{{))+Pe{l)]/de, respectively. 

For each disorder configuration, we can calculate the 
energy levels £j from the scattering matrix S{£,x') of a 
wire with Hamiltonian ([T]) for x < x' and H = v^pa^ 
for X > x'. The relation between S (which is really a 
complex number of unit modulus in the present case) 
and the energy levels is given by the Friedel sum rule. 



(a) 



N{e)^{2TTir^ lim \ndet[S{e,x')S{e,-x'y 



(3) 



where IuUx'^-oo S{e,x') = i. The probability distribu- 
tions Peif) can then be obtained by considering the evo- 
lution of S{e^x') upon repeatedly increasing x' by the 
small amount 5x' <C min(Z,^), ^ = tip/A being the super- 
conducting coherence length in the topological domain. 

The evolution of 5(£, x') takes the form of a Langevin 
process This Langevin process takes its simplest 

form if we use the parametrization 



S — i tanh y, 



(4) 



instead of the standard parametrization of S through the 
scattering phase, as the parametrization Q makes the 
noise term become independent of S. The parameter 
y takes values on the real axis ±i7r/4, see Fig. [T^, and 
is continuous at y = ±oo. Concatenating the scatter- 
ing matrix for a given x' with the scattering matrix of an 
added slice of length Sx' , we can compute the correspond- 
ing change Sy of y. This yields the Langevin process 

{Sy) = — [zssmh2y-A{x')], {5y^) = °^. (5) 

Near y — ±oo, the shifts are dominated by the term pro- 
portional to £, which unidirectionally couples the branch 
at y = ±oo=Fi7r/4 into the branch at y = ±oo±i7r/4, see 
Fig. [l] The initial condition for x' < is y = oo -I- Z7r/4. 
For x' > L the Langevin process returns y to the starting 
point y = oo-|-i7r/4. The return takes place via the upper 
branch at lm.y = -k /A if Imy(L) = 7r/4, and via the lower 
branch otherwise. 

The Friedel sum rule ^ now identifies Pe{N) as the 
probability distribution of the number of times N that 
the variable y has passed through the point at — oo upon 
increasing x' from a;' = to a;' = L. We have calculated 
this probability distribution through direct numerical 
simulation of the Langevin process, as well as through an 
asymptotic analysis valid in the limit s <Si min(A, l/r). 

For the asymptotic analysis, we observe that for e ^ 
A the Langevin process ^ is dominated by the term 
proportional to e if |Rei/| > yi^), with 



1 max(2A, l/r) 
= 2 1^ 



(6) 
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FIG. 1: (color online) (a) The variable y in the parametriza- 
tion Q of the scattering matrix takes values on the real axis 
±in/4. The black dotted arrows indicate the boundary con- 
ditions at y = ±oo. In the simplified Langevin process used 
for the asymptotic analysis ([7|, the boundary conditions are 
at Rey = ±y (full black arrows). The direction of the drift 
term in topological {0 < x' < L) and non-topological {x' < 
and x' > L ) regions is indicated by the full red arrow and 
dashed red arrow, respectively, (b) Typical dependence of the 
Majorana end-state energy eo and the energies ei, £2, ■ . . , of 
low-lying bulk states on the system size L. 



whereas the other terms dominate if |Rey| < y(e). To 
logarithmic accuracy, we approximate the Langevin pro- 
cess (13) by a truncated process in which |Rey| < y 
and the energy contribution to {6y) is omitted 24J. The 



unidirectional connection between the upper and lower 
branches now takes place at Rey = ±y{e), see Fig. [I|d. 
The resulting Langevin process for y is then specified by 
the equations 



{6y) = -5x'/i, {Sy^} = 5x' /2l 



(7) 



with absorbing boundary conditions ("sink") at y = 
±(y(e) — j7r/4) and hard wall boundary conditions at 
y = ±(y(e) + W4), see Fig. [TJd. 

There is a qualitative difference between the Langevin 
processes at the upper branch (Imy ~ 7r/4) and the lower 
branch (Imy = — 7r/4): At the upper branch, the drift 
term proportional to 1/^ pushes the variable y towards 
the sink, whereas at the lower branch it keeps it away 
from the sink. The slow diffusion in the latter case does 
not affect Eq, but it dominates the probability distribu- 
tion of all higher levels Sj, j = 1, 2, . . .. By analyzing the 
diffusion process on the upper branch, wc find that the 
probability P^iO) is [25J 



PeiO) 



1 



erfc 



L/^ - 2y(£) 



(8) 



Recalling that y{s) is given by Eq. (|6]), we conclude that 
ln(£o/2A) has a normal distribution with mean and vari- 
ance given by 



(ln(£o/2A)> = -L/t 
var ln(eo/2A) = L//, 



(9) 



up to corrections of order unity that can not be de- 
termined from the above argument. Similarly, analyz- 
ing the diffusion process on the lower branch, we find 
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Pe{0) + Ps{^) = e ^^('^)), where the disorder-averaged 
integrated density of states {N{e)) cx (L/0(e/A)2'/« 
[20l El] , with proportionahty constant that could not be 
determined from the asymptotic analysis. From this, we 
find that pi(e) oc (Lr/0(e/A)2V?-i for small e. The 
theoretical predictions are compared to numerical simu- 
lations of the Langevin process in Fig. [2] 

One- dimensional spinless p-wave superconductor. — We 
now extend these results to a one-dimensional p-wave su- 
perconductor, using a continuum version of the model 
considered in Refs. HI [51 



phase (j) decouple, 

6x' Sx' 
{Sy±) = — (iesinh2y± - A) + — coth(2/± - yzp), 
vp 21 

{5yl) - -{5y±Sy^) = 6x'/2l, 

{54>) = 2pfSx', 

{S<j>^) = A5x' /l + {5L/2l)coih^{y+-y-). (13) 



The initial condition is 



H = 



P 

2m 



+ V{x) - f?j a- 2 - A'pcr^. 



(/)(0)=0 and y± (0) = ±cx3 ± i7r/4. 



(14) 



(10) 



For a given disorder realization, the solutions Sj of Eq. 



Here /x is the chemical potential, V{x) the disorder po- 
tential, which we take according to the same distribution 
as in the previous case, and cr^ and are Pauli matrices 
acting in the electron-hole space. We model a topological 
domain of finite length L by setting /i = —oo for x < 
and X > L, and n = p^ /2m > for < x < L, pp — mwp 
being the Fermi momentum |26) . Further, A = A'pp is 
the superconducting gap. Throughout our calculation we 
assume that fi/pF- 



(12) oscillate as a function of L, with oscillation period 
~ tt/pf, as shown schematically in Fig. This follows 
from the observation that and y- are "slow" as a 
function of L [they vary on the scale min(^,^)], whereas 
(j) is a "fast" variable {6(j)/6x' « 2pp): Starting from the 



initial condition (14), solutions of Eq. (12) then appear 



Previous studies of lattice versions of the model ( 10 ) 
addressed the disorder-averaged density of states (i^(£)) 
in the limit L — >■ oo [SJiMj. Using a strong-disorder renor- 
malization group approach, Motrunich et al. showed that 



in quick succession upon increasing L at fixed e, until y-^- 
passes through the point at — oo or y_ passes through 
the point at -t-oo, whichever occurs first. No solutions 
of Eq. ( 12 ) are found upon increasing L further, until 



eventually again one of the variables y± passes through 



the model ( 10 ) is in a topological phase if the disorder g^^^ . ^j? ^-^g 



a point at infinity and solutions to Eq. (12) reappear in 
quick succession, cp. Fig. [TJd. 

We now calculate the probability distributions po.max 



strength is below a critical value, and in a non-topological 
phase for stronger disorder. On both sides of the critical 
disorder strength, the density of states J/(e) has a power 
law dependence on e for e ^ A, with an exponent that 
depends on the disorder strength. For the continuum 
model ( 10 ) we now show that the transition is at ^ = 21 



maximum £o,i 



and the minimum 
£i,min with respect to variations of L of order n/p-p. (Note 
that EOjinax and ei^min are the energies relevant for setting 
the operation speed of a hypothetical topological quan- 
tum computer.) Repeating the arguments of the first 
part of this paper, these probabilities obey 



and calculate the probability densities of the Majorana 
end-state energy eg and the lowest-lying bulk level for 
disorder strengths below the critical value. 

Our calculation essentially follows the approach taken 
above for the Dirac equation with random mass with 
some modifications. We define a 2 x 2 scattering ma- 
trix 5(6, x') of a wire with Hamiltonian given by Eq. 
(10) for X < x' and by H — {p^/2m)az for x > x' and 



d 



]'<2j 



Pj,m\n 



de 



(15) 



(16) 



J'<2j-1 



parametrize S through 



S{e,x') 



2^. 



±e*'^tanhj/± — itanhy± 
ztanhj/± ±e~"^tanhy± 



(11) 

where the variables and y_ take values on the real axis 
±i7r/4, see Fig.[lJ and is a real phase. The energy levels 
can no longer be calculated from the Friedel sum rule 
([3]), but instead have to be obtained from the condition 
dct[l + S{e, L)] = 0, which becomes 



cos — coth(y_ — j/-|_). 



(12) 



in the parametrization (11). For h/pp 5x' <^ /,^, the 



resulting Langevin processes for the variables y± and the 



where Pe{N) is the probability that (in total) the vari- 
ables y+ or y_ have passed N times through the points 
at ±00 upon increasing x' from to L. We have calcu- 
lated these probabilities from direct numerical simulation 
of the Langevin process, as well from an asymptotic an- 
alytical solution valid in the limit e ^ min(A, 1/t). For 
the asymptotic analysis, we make the same simplification 
of the Langevin process as in the case of the Dirac equa- 
tion with random mass. In addition, we observe that 
for the energies of interest, one of the variables y± effec- 
tively remains pinned at —y{s) — in/A, so that the factor 
coth(j/+ ^2/-) in the interaction term may be approx- 
imated by ±1. The resulting Langevin process for the 
remaining variable is then specified by the equations 

{Sy)^-SL/^ + SL/2l, {5y^) ^ 5L/21, (17) 



4 
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small in L/^. This implies that in principle, there is a 
large parameter window in which both conditions on the 
operation speed can be met if L is made sufficiently large. 
It is important to note, however, that with increasing dis- 
order or increasing L, this parameter window is shifted to 
lower energies which would require the topological quan- 
tum computer to operate at a lower temperature and 
lower speed. 
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and the Alexander von Humboldt Foundation. 



FIG. 2: (color online) Integrated probability density Pe(0) 
of the Majorana end-state energies £o (blue dots) and £o,max 
(red crosses) in the two models obtained from a numerical 
solution of the Langevin process, together with the theoretical 
prediction (solid). Energies are normalized to the median em 
of the distribution. Left inset: Logarithm of the integrated 
probability distribution Ps{Q) + fe(l) of the lowest lying bulk 
state energies e\ and ei.min. Right inset: Average of logeo 
and logeo, max vs. length. In both insets the slope of the solid 
lines is given by the theoretical predictions of the main text. 



with the boundary conditions as specified below Eq. ([t]) . 
The result for Pe(0) has the same functional form as in 
the case of the random-mass Dirac equation, and we con- 
clude that ln(eo,max/2A) has a normal distribution with 
mean and variance given by 



(ln(eo,n: 
var ln(eo^ 



c/2A)> 
.x/2A) 



-L[i/e- 

L/21, 



1/(20], 



(18) 



up to corrections of order unity that can not be deter- 
mined from the asymptotic argument. The end-state en- 
ergy remains exponentially small in L as long as 21 > ^, 
which identifies 2Z = ^ as the critical disorder strength 
that drives the system into the non-topological phase. 
Similarly, we find Pe(0)-FPe(l) = e-<^('^)> , with {N{e)) cx 
(L/C)(£/A)4V?-2 [51|Ttl[2g, from which we conclude that 
Pi,min(£) oc (Lr/^)(£/A)'*'/^~'^ for small energies. At the 
critical disorder strength, the integrated density of states 
takes the Dyson form N{e, L) cx ln^(£/A) [SinillT]. The 
theoretical prediction for po,max is compared to numerical 
simulations of the Langevin process in Fig. [2j 

Conclusions. — For both models of a topological super- 
conducting wire, we find that the energy splitting eo of 
the Majorana end states has a log-normal distribution, 
implying large sample-to-sample fluctuations. Neverthe- 
less, for sufficiently long wires, the width of the log- 
normal distribution remains small compared to its av- 
erage. In this case, the lower limit on the speed of the 
qubit operations is well determined by the typical value of 
the log- normal distributions in Eqs. ^ and ( 18 ), which is 
exponentially small in L/^. By contrast, we find that the 
energy ei of the lowest-lying bulk state is algebraically 
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